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ABSTRACT 

Eclipsing binary DI Herculis (DI Her) is known to exhibit anomalously slow apsidal precession, below 
the rate predicted by the general relativity. Recent measurements of the Rossiter-McLauglin effect 
indicate that stellar spins in DI Her are almost orthogonal to the orbital angular momentum, which 
explains the anomalous precession in agreement with the earlier theoretical suggestion by Shakura. 
However, these measurements yield only the projections of the spin-orbit angles onto the sky plane, 
leaving the spin projection onto our line of sight unconstrained. Here we describe a method of 
determining the full three-dimensional spin orientation of the binary components relying on the use 
of the gravity darkening effect, which is significant for the rapidly rotating stars in DI Her. Gravity 
darkening gives rise to nonuniform brightness distribution over the stellar surface, the pattern of which 
depends on the stellar spin orientation. Using archival photometric data obtained during multiple 
eclipses spread over several decades we are able to constrain the unknown spin angles in DI Her with 
this method, finding that spin axes of both stars lie close to the plane of the sky. Our procedure fully 
accounts for the precession of stellar spins over the long time span of observations. 

Subject headings: (stars:)binaries:eclipsing — stars:rotation 



1. INTRODUCTION. 

The binary system DI Herculis (DI Her, HD 175227) 
was discovered as an eclipsing variable by Hoffmeister 
(1930). It consists of two massive B stars on an eccen- 
tric orbit (e = 0.49) with period P = 10 d .55, inclined 
at an angle i = 89.3° with respect to our line of sight 
(see Table Q] for parameters of both binary components). 
A unique property of this system that has been attract- 
ing a lot of attention for almost three decades is its low 
rate of apsidal precession, cj f, s = 1.24° ± 0.18°/100 yr 
(Martynov & Khaliullin 1980). This is almost two times 
lower than the general relativistic apsidal precession rate 
wgr = 2.35°/100 yr theoretically predicted based on the 
measured parameters of the system (Rudkjobing 1959). 
For a long time this discrepancy was not understood 
(Maloney et al. 1989) and even ascribed to the failure 
of general relativity (Moffat 1989). 

Shakura (1985) suggested that slow apsidal precession 
in DI Her is caused by the misalignment between the 
spin and orbital angular momentum axes of the sys- 
tem. Indeed, for stellar spins strongly misaligned with 
the orbital angular momentum the rotation-induced stel- 
lar quadrupole gives rise to a contribution to to with a 
sign opposite to that of wgr- Then, if this quadrupole- 
induced precession is fast enough it can easily alter the 
full rate of apsidal precession and even make it smaller 
than wgr as in the case of DI Her. Somewhat less ex- 
treme version of this idea has recently been applied to 
another eclipsing binary AS Camelopardalis (Pavlovski 
et al. 2011), which also exhibits relatively slow rate of 
apsidal precession. 

The spin-orbit misalignment in DI Her has been re- 
cently confirmed by Albrecht et al. (2009) who used 
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the evolution of stellar spectral signatures during the 
eclipse (the so-called Rossiter-McLaughlin effect, Holt 
1893; Rossiter 1924; McLaughlin 1924) to set constraints 
on the spin orientation of both stars. They found that 
both stars of DI Her have their spin axes nearly perpen- 
dicular to the orbital angular momentum axis, which is 
at odds with the common wisdom regarding spin-orbit 
orientation in close binary stars, but can naturally ex- 
plain the slow apsidal precession in this system. 

Unfortunately, the Rossiter-McLaughlin effect allows 
one to determine only the sky plane projection A of the 
angle a between the spin and orbital momentum axes for 
the stars. The angle j3 between the stellar spin and our 
line of sight remains essentially unconstrained. However, 
figuring out whether spin-orbit misalignment can explain 
the observed ui does depend on the value of (3, since the 
stellar spin frequency uj is inferred from the measured 
projected stellar rotation speed tVotsin/3 = wi?*sin/3 
(R+ is the stellar radius). Albrecht et al. (2009) and 
Claret et al. (2010) tackled the issue of undetermined 
j3 by means of Monte Carlo simulations, assuming this 
angle to be uniformly distributed. 

In this work we develop a method of analyzing the pho- 
tometric eclipse data, which allows us to constrain the 
angle /3 without using spectroscopic data. This method 
relies on the fact that both components of DI Her are 
rapidly rotating stars (iVot sin /? exceeds 100 km s _1 for 
both components, see Table []}, and must exhibit a non- 
uniform surface brightness distribution due to the gravity 
darkening effect (von Zeipel 1924). This surface bright- 
ness pattern is sensitive to the orientation of stellar spin 
axis with respect to our line of sight. By probing the 
brightness distribution via the detailed shape of the sys- 
tem lightcurve during the eclipse one can infer the full 
spin orientation of both stars. Analogous method was 
recently proposed by Barnes (2009) for analyzing plane- 
tary transits around rapidly rotating stars, and applied 
by Szabo et al. (2011) and Barnes et al. (2011) to deter- 
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Fig. 1. — Observer coordinate frame on the left, and a sketch of 
the orbital geometry on the right. See text for more details. 

mine the spin-orbit misalignment in a transiting system 
KOI-13.01. 

In determining the unknown angles P for both stellar 
components we also use other constraints on the system 
parameters, such as the observed apsidal precession rate 
and the evolution of the projected stellar rotation veloc- 
ities v ro t sin p over the long time span. 

This work is organized as follows. In $2] we describe 
geometric setup of the problem, and eclipse lightcurve 
modeling. In fJ3]we describe observational data and our 
fitting procedure. Our results are presented in g] and 
discussed in fj5] In Appendix [B] wc derive equations de- 
scribing evolution of the system orientation as a result of 
spin and orbital precession, which may find other appli- 
cations. 

2. ECLIPSE MODELING. 
2.1. Geometry of the system. 

To model eclipse lightcurve we use two Cartesian co- 
ordinate systems. One is the observer frame (X,Y,Z), 
which describes the orbital orientation of the system: Z 
axis points from the system barycenter towards the ob- 
server, X axis is along the direction of the sky-projected 
orbital angular momentum and Y axis is aligned with 
the line of nodes, see Figured] 

Another system (xi, Ui,Zi), where i = p,s for primary 
and secondary, respectively, is aligned with the stellar 
symmetry axis {symmetry frame): its Z\ axis is along the 
stellar spin angular velocity ujj, and Xi and yi axes are 
obtained from X and Y by two rotations: first, a rotation 
around Z axis by the angle Xi and then another rotation 
around axis obtained from Y in previous step by the 
angle ft. We will use the symmetry frame to describe 
the stellar surface shape and temperature distribution 
and the observer frame to characterize the visible sky- 
projected stellar disc. 

Spin angular velocity <jjj in the observer frame is then 
given by 

(sin fa cos A A 
sin Pj sin Aj (1) 
cos Pj J 

The vector towards the observer in the symmetry sys- 
tem of each star is 



(2) 



To simulate the eclipse light curves we also need to 




know the relative stellar trajectory projected onto the 
plane of the sky. The projected position of the center 
of the secondary with respect to primary in the observer 
frame is given by 



(— cos«sin(/ + vj) 
_ cos(/ + w) 
sini sin(/ + zu) 



(3) 



where r s is the distance of the secondary star with respect 
to the primary and / is the true anomaly (see Murray & 
Dermott 2000 for relation with other orbital parameters) . 
Equation Q fully determine the time evolution of R s 
during the eclipse. 

2.2. Shape of the stellar surface. 

We now describe the variation of intensity of emission 
over the stellar surface due to the gravity darkening ef- 
fect. Our results will be valid both for the primary and 
the secondary, so we will omit the subscript i = p, s in 
this subsection. 

First, we determine the geometry of the sky-projected 
disc by assuming the shape of the stellar surface to coin- 
cide with isopotential surfaces for the effective potential 

4> eff (x, y,z) = - GM : - \tll{x 2 + y 2 ), (4) 
yr + y^ + z z L 

as a function of coordinates in the "symmetry" frame. 
Assuming slow rotation one can obtain the equation for 
the stellar surface in the form 



x2 + y 2 +z i = B 2 



(5) 



where R po i is the polar radius of the star. Thus, because 
of the rotation stellar surface has an ellipsoidal shape 
with oblateness r\ given by 
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(6) 



Here the parameter S is related to the ratio of the 
rotation rate u)i to the breakup rotation rate uj^ = 
(GM-,,/ R^) 1 / 2 at which the centrifugal force balances 
gravity at the stellar surface. Equation (JS|) can be 
re- written in polar coordinates as r(6) = R po i e (l + 
S sin 2 6»/2). 

The actual angular velocity of stellar spin is calculated 
from the spectroscopically measured projected stellar ro- 
tation velocity (v ro t sin P) b s as 



(vrot smP) obs 
R eq sin/3 



(7) 



where we assumed R eq equal to the stellar radius quoted 
in the literature. 

2.3. Temperature distribution over the stellar surface. 

Because of rapid stellar rotation the brightness tem- 
perature of the stellar surface is not constant but obeys 
the von Zeipel (1924) law 



T(R) — T po i 



g C ff(R.) 

9 P oi 



Tpoiip(B.), 



(8) 
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Fig. 2. — Latitudinal distribution of the brightness temperature 
for a stellar model with rotation parameter S = 0.09 (points; De- 
upree, 2010) fitted by a von Zeipel law (equation JS}; solid line) 
with low value of (3 g = 0.075 that is close to /3 g = 0.1 used in this 
work. 

where T po i is the value of T at the stellar pole, (3 g is 
the gravity darkening power law index, R is the three- 
dimensional radius-vector from the stellar center to a 
point on the stellar surface, and g c g is the local effec- 
tive gravitational acceleration: 

GM 

gefl(R) = — + oj 2 R ± . (9) 

Here R^ = R — (u> ■ R) /to is the distance to that point 
from the stellar spin axis, and g po i is the value of g e ff at 
the stellar pole, where R± = 0. 

Conventional gravity darkening theory (von Zeipel 
1924) predicts f3 g — 0.25. However, recent detailed theo- 
retical calculations (Deupree 2011) of the latitudinal dis- 
tribution of the effective temperature for rotating stars 
performed in the wide range of stellar masses (between 
1.625 and 8 M ) suggest a considerably weaker depen- 
dence of T(R) on g c g. We illustrate this point in Figure 
[H where we plot the latitudinal distribution of the ef- 
fective temperature for a particular stellar model from 
Deupree (2011) corresponding to a rotation parameter 
£ = 0.09. This distribution depends on stellar mass only 
weakly, meaning that it can be applied for DI Her com- 
ponents as well. One can see that equation (|8]) with 
fig = 0.075, which is considerably lower than 0.25, pro- 
vides excellent fit to these data. 

On the observational side, interferometric measure- 
ments for rapidly rotating stars by Che et al. (2011) 
find P g ss 0.146 for 1.77M© star /3 Cassiopeiae, having 
v ro tsm/3 ss 75 km s _1 and (3 g « 0.19 for 4.15M Q star a 
Leo, rotating at iv ot sin/3 « 340 km s -1 . Even though 
the latter is very similar in mass to the DI Her compo- 
nents, it spins much faster (spin parameter S is almost an 
order of magnitude higher), making direct extrapolation 
to the DI Her case difficult. Despite these ambiguities, it 
is clear that both theoretically and observationally one 
typically infers fid < 0.25. In this work we have chosen 
to adopt jig = 0.1 more in line with the work of Deupree 
(2011). 

2.4. Intensity distribution over the sky-projected disc. 

Equations ©, © provide us with a simple expres- 
sion for T(R) in the symmetry frame of a star, since 
R in this frame can be trivially derived from equation 
([5]). However, for the purposes of eclipse lightcurve 




Fig. 3. — Intensity distribution for the primary (a) and secondary 
(c) stars, and temperature distribution for the primary (b) and 
secondary (d). Calculations assume j3 v = 70°, \ p = 72°, /3 S = 
110°, A s = -84°, S = 0.040 for the primary star, S = 0.046 for 
the secondary star, limb-darkening coefficient ci = 0.35 for the 
primary star, c\ = 0.64 for the secondary star and von Zeipel 
parameter /3 g = 0.1. The dot on each plot indicates the position of 
the stellar pole. The blue line shows the trajectory of the center of 
the secondary (primary) star during primary (secondary) eclipse. 
The arrow shows the direction of the projected orbital angular 
momentum L. Temperature distribution clearly illustrates gravity 
darkening. The intensity variation is large (tens of per cent) mainly 
because of the limb darkening effect, which is much larger than the 
gravity darkening. 

modeling we need to know the temperature distribu- 
tion in the observer frame, projected onto the plane of 
the sky. In Appendix [S] we describe the relation be- 
tween (xi,yi,Zi) and (X,Y,Z) frames, which allows us 
to write the dimensionlcss function ip in equation (|8|) as 
V>(R) = i/>(X, Y, Z(X, Y, o>i)) = ip(X, Y, uj % ) for each star. 
The dependence of ip on u>j is the key factor that allows 
us to use stellar photometry during eclipse to determine 
the spin orientation of the stars. Figure \3jp,d illustrates 
the distribution of the brightness temperature over the 
stellar surface projected onto the sky plane. 

Spectral density of the stellar radiation flux detected 
on Earth is 

Fx = h J h{T P ort{X,Y,u>))® x {X,Y,u;)dXdY, (10) 

5 

where d is the distance to the system, $a is the limb- 
darkening law, and I\(T) is the spectral intensity at 
a given temperature T. In this work we assume that 
I X (T) = B\(T), where B\(T) is a standard black-body 
radiation function. Thus, spectrum of each star is in gen- 
eral a multi-color blackbody, parametrized by T po i and 

On the other hand, values of the effective temperature 
for both stars quoted in the literature are obtained as- 
suming that both stars radiate as pure blackbodies char- 
acterized by a single value of temperature, uniform across 
the stellar surface (Claret 2010). In that case the total 
bolometric flux is F = L»/(47rcP) = aT^^R^/d 2 , where 
L» is the stellar luminosity. This assumption is not valid 
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for rapidly rotating stars because of gravity darkening, 
and T po i cannot be taken equal to T*. To relate them 
we integrate equation (fT0|) over all wavelengths to obtain 
the following expression for T po i\ 



T, 



pol 



^ / ip A (X,Y,u>)$ x (X,Y,u>)dXdY 



-1/4 



(11) 

i.e. there is a correction depending on the orientation of 
stellar spin. This self-consistent derivation of T po i is an 
important part of our procedure which distinguishes it 
from the approach of Barnes (2009). 

For simplicity the limb-darkening law in this work is 
assumed to be frequency- and spin-independent and have 
a simple functional form 



$a = 1 - ci(l - fi), 



(12) 



where c\ is constant and /i is the cosine of the angle be- 
tween the local normal £ to the stellar surface and the 
observer's line of sight: fj, = £ ■ n. The total measured 
flux in the V band Fy is obtained by additionally con- 
volving F\ in equation (|10[) with W\ — the normalized 
transmission function for that band — over A. 

Figure [3K,c shows how the radiation intensity is dis- 
tributed over the stellar surface for the two components 
of the DI Her system out of eclipse, when both gravity- 
darkening and limb-darkening are taken into account. It 
is obvious that limb-darkening has a much stronger effect 
on the intensity distribution than the gravity darkening, 
complicating measurement of the latter effect in the pho- 
tometric data and determination of the spin orientation 
of the two stars. On the other hand, as long as the stel- 
lar spin axis is not aligned with our line of sight, the 
gravity darkening results in a non-axisymmetric bright- 
ness distribution with respect to our line of sight, unlike 
the limb darkening, which is axisymmetric. This helps 
one disentangle the two contributions in the photometric 
data. 

For simulating the light curves we need to calculate the 
flux blocked during the eclipse 



/,.,;. , (w. t) = J H(X, Y, t)F v<i (X, Y, u)dXdY, i = p, s, 

(13) 

where H(X, Y, t) equals 1, if the secondary (primary) star 
blocks starlight of the primary (secondary) at position 
(X, Y), and if not. Then the total flux observed on 
Earth is 

/(a>,t)=Jp+J a -/ M , i (a»,t), (14) 

where I p and / s are the unblocked stellar fluxes, i.e. cal- 
culated from equation (|13|) with H set to unity. During 
the eclipse Ibi t i(u},t) changes because function H varies 
with time across the surface of eclipsed star. In this work 
we take the surface (and frequency because of the finite 
bandwidth) integral in equation (fl"3)) using Monte Carlo 
technique. 

Our lightcurve modeling procedure is illustrated in Fig- 
ure IU where in left panels we demonstrate the differences 
between the model assuming uniform distribution of the 
surface temperature and the models in which the effect 
of gravity darkening is fully taken into account. One 
can see that for relatively fast rotation (S = 0.13 or 
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Fig. 4. — (Left) Top panel shows simulated light curves of the 
secondary eclipse for the uniform temperature distribution (solid 
line), fast rotation (S=0.13) and two different secondary spin ori- 
entations - /3 s i = 70°, A s i = 90° (dashed line) and s2 = 70°, 
A s 2 = —100° (dotted line) (primary star parameters are kept con- 
stant at f3 p = 70°, A p = 72°). On the lower left panel the difference 
between simulated curves (dashed for the first case and dotted for 
the second one) and the case of uniform temperature distribution 
is given. For the first orientation dimmer regions are blocked first, 
with the opposite being true for the second one. (Right) Top panel 
shows simulated light curves for the uniform temperature distribu- 
tion (dashed line) and actual DI Hercules parameters (/3 P = 70° , 
0s = 110°, A p = 72°, A 3 = -84°, S p = 0.040, S s = 0.046). Dots 
correspond to observational data (secondary eclipse 7/13/1986). 
The differences between the theoretical lightcurves are hardly visi- 
ble, so we visualize them on the lower panel (we do not show data 
point there as they would be off scale). 

uj = 0.36tdb) the difference between the uniform T case 
and the gravity-darkened models is at the level of sev- 
eral per cent, which should be easily detectable in single- 
epoch observations. 

In right panels of the same Figure we show the compar- 
ison between the uniform temperature case, the gravity- 
darkened model with spin angles resulting from our fits 
to the data (see In this case S « 0.04, see Table 
[H and the difference between the uniform and gravity- 
darkened models is very small, ~ 10~ 3 mag. Thus, in the 
case of DI Her one would need very high-quality photom- 
etry (at the level of several x 10~ 4 mag) to detect gravity 
darkening-related asymmetries in the lightcurve shape. 

3. OBSERVATIONAL DATA AND FITTING PROCEDURE 

3.1. Observations. 

The dataset we use for determining spin orientation of 
the DI Her components consists of photoelectric V band 
measurements of this system using the 50-cm AZT-14 re- 
flector at the Tien Shan Observatory of the Astrophys- 
ical Institute of Kazakhstan and the Zeiss-600 reflector 
at the Crimean Station of the Sternberg Astronomical 
Institute prepared in 20032008. The database also con- 
tains the photoelectric observations going further in the 
past by Semeniuk (1968), the 1968-1978 observations by 
Martynov and Khaliullin (1980), the 1986-1988 observa- 
tions by Khodykin, Volkov, and Metlov, and the 2004 
observations by Shugarov (see Kozyreva & Bagaev 2009 
and the references there in). These observations use dif- 
ferent instruments, which were not cross-calibrated. Our 
fitting uses 9 eclipses, 4 primary and 5 secondary. The 
photometric errors are unconstrained in all cases and we 
describe in $3A\ how we deal with this issue. 



3.2. 



Time evolution of the system. 



High mass, rapid rotation and relative proximity of 
the stars in DI Her system drive rapid evolution of the 
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Fig. 5. — Eclipse data used in our modeling with best fits overplotted. Horizontal axis is JD and the vertical one is apparent magnitude 
in the V band for every eclipse. Labels indicate the date of eclipses and the adopted noise levels. 

cates the eclipse fitting, but on the other hand it provides 
us with a unique opportunity to use the signatures of this 
variation of stellar spins in eclipse modeling on long time 
intervals. 

In this regard our procedure of eclipse simulation uses 
somewhat different strategy that the one proposed by 
Barnes (2009): instead of using high-accuracy photomet- 
ric data in a single epoch, which provides sensitivity to 
spin orientation only through the shape of the eclipse 
lightcurve, we use low-accuracy photometry obtained at 
different epochs, which allows for an additional effect of 
the current spin orientation of DI Her on the fitting pro- 
cedure — through the time evolution of cjj. 
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Fig. 6. — Evolution of the DI Her spin orientation (angles f) 
and A) for both stars and the inclination i of the system both on 
short (left) and long (right) time intervals. Dashed vertical lines 
correspond mark the locations of eclipses used in our modeling. 
Variation of spin angles ft and A;, i = p, s are very significant. 

spin orientations for both components. In Appendix [B] 
we derive equations that describe evolution of the stellar 
spins and orbital elements of the system. In particular, 
we show there that stellar spins in DI Her can rotat e 
by more than 100° within a century, see equation (IB14[) . 
This evolution obtained by integrating equations from 
Appendix [Bl over time is illustrated in Figure |H1 

This Figure clearly demonstrates that spin orientation 
of the DI Her components significantly changes over the 
time span of our full dataset. On one hand this compli- 



3.3. Evolution of v rot sin /3 . 

Precession of stellar spins causes evolution of 
v rot _i sin Pi (i — p, s) on long time interval, which can 
be compared with past measurements. We took the 
values of u ro t,i sin /3j from Albrecht et al. (2009), who 
compiled measurements from different epochs starting in 
1948. Figure ([TD]) shows these data. Unfortunately, the 
large (or completely undetermined as in the case of 1948 
data point) uncertainties of these measurements do not 
allow us to infer the values of j3 p and /3 S from these data 
alone. However, even the weak constraints based on these 
measurements turn out being quite useful. 

We will use the fact that based on these data v rot sin /3 
is currently increasing for both primary and secondary 
components. As v rot is constant this can only be due to 
sin j3j increasing in time. From the evolution equation 
(|B20I) we find the expression for the derivative of sin /3j 
U = P, s) 



dsm-Pj . . . 
; — - = 61 1 p i sin i sin A, cos 

dt 3 3 

3 

-fipj sin 2 i sin(2Aj) cos 



dt 

- sin i sin /3j cos Xj ) ~ 



■ j3j (cos i cos f3j 



i + 

i&,(15) 
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where flp.j is the frequency describing spin precession 
caused by the rotation-induced stellar oblateness (see 
equation (|B14I0 and the approximation holds for i m 
90° . From the Rossiter-McLaughlin measurements of Al- 
brccht et al. (2009) we know that sin 2Xj is positive for 
the primary and negative for the secondary. Thus, it fol- 
lows from the current time derivatives of sin f3 PtS that /3 P 
must be less than 90° while /3 a should be greater than 
90°. We use this information to help constrain stellar 
spin orientation in jQ 

3.4. Photometric fitting procedure. 

The specific parameters of DI Herculis that we used in 
simulations are summarized in Table [TJ To keep things 
simple we have only varied the most important unknown 
quantities — the angles /3 p and /3 S . Given the quality of 
the data we expect that fitting for extra parameters in 
our model would result in too many degeneracies between 
the different model variables. 

We integrate back in time the evolution equations for 
spin and orbital parameters (Appendix |Bj starting from 
15 July, 2008, which is set as the initial point in our cal- 
culations. The angles j3 p and j3 s that we vary correspond 
to this epoch. The other two angles specifying spin ori- 
entation were fixed at X p — 72° and X s = —84° based 
on the Rossiter-McLaughlin measurements of Albrecht 
et al. (2011). 

To obtain better eclipse fits we had to introduce quite 
different limb-darkening coefficients for the two compo- 
nents, which, given the proximity of stellar masses in DI 
Her, suggests that our data are affected by some system- 
atic effects. Nevertheless, given that the limb darkening 
only weakly affect the non-axisymmetric surface bright- 
ness distribution due to gravity darkening (see ^2.41) , the 
actual values of limb-darkening coefficients are not so im- 
portant. 

We constrain DI Her spin orientation as follows. For 
each pair /3 P , j3 s we compute 
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(16) 



where index j runs through all 9 eclipses, i runs through 
the number of data points per each lightcurve (total of 
Nj for j-th eclipse), N = Y2j=i IS the total number of 

data points, I(t) is given by equation (fl"4|) . I 3 obs {t) is the 
observed intensity, and o~i is the variance. The best fit 
values of (3 P and j3 s are determined by finding the mini- 
mum of x 2 over a large two-dimensional grid of values of 
these angles. We use only the eclipses with well-defined 
minima. We did not try to match theoretical and obser- 
vational eclipse minima with our direct backward inte- 
gration in time and instead just shift theoretical curves 
horizontally by small amount at each epoch for a better 
fit. 

As mentioned before, the errorbars for our dataset are 
not constrained, so we employed the following procedure 
to estimate them. First, we took all cr, to be constant 
and run our minimization procedure to find the best fit 
values of (3 p and j3 s - Second, for each out of 9 eclipses we 
measure the scatter crj of the observational data points 
around the model lightcurve computed assuming these 
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Fig. 7. — Sensitivity of spin angle determination to the noise level 
assumed for the data. We sim ulat e the primary eclipse lightcurve 
using theoretical prescription H14H with a given level of the Gaus- 
sian noise a for a fixed spin angle /3 P = 40° for the primary. We 
vary only this angle in our \ 2 fitting (/3 S = 40°, A p = 72° and 
\ 3 = —84° are constant here) just to illustrate that large noise level 
does not allow us to constrain system parameters, while a a, 0.01 
mag yields the correct orientation of the system. 

particular values of /3 p and /3 S . This provides us with 9 
different values of sjj (indicated in panels in Figure [5] for 
each eclipse), which we use as error estimates in equation 
p^|) . Typical values of o~j found using this procedure are 
~ 0.01 mag. We then perform the final \ 2 minimiza- 
tion adopting these values of o~j as error estimates for 
corresponding eclipses. In this approach all data points 
corresponding to j'-th eclipse have a single value of the 
photometric error equal to o~j. 

We test the performance of our fitting algorithm by ap- 
plying it to a simulated dataset in a simplified setup. We 
calculate a theoretical primary eclipse lightcurve includ- 
ing the gravity darkening effect and assuming a binary 
with physical parameters (M*, R*, T*, etc.) of the DI 
Her (in particular with S ~ 0.04 for both stars). We 
take somewhat arbitrarily (3 p = 40°, j3 s = 40°, \ p = 72° 
and \ s = —84°. We then add some random Gaussian 
noise with variance a to this simulated lightcurve. For 
this test we assume /3 S , X p and A s to be known and try 
to measure the value of only /3 P using our procedure. As 
a consequence, we need to perform only one-dimensional 
minimization over f3 p . 

The results of this exercise are shown in Figure [7) 
where we show \ 2 curves for two different values of the 
noise variance a: 0.01 mag and 0.03 mag. One can see 
that for st = 0.03 mag our parameter estimation proce- 
dure cannot recover the adopted value of f3 p — the x 2 dis- 
tribution has very extended flat bottom which does not 
lead to a useful constraint on j3 p . However, for st = 0.01 
mag our procedure works reasonably well and the min- 
imim of x 2 is close to the input value of j3 p — 40° . Since 
the simulated lightcurve was computed for realistic phys- 
ical parameters of the DI Her and the noise levels for in- 
dividual eclipses in real data o~j are indeed ~ 0.01 mag we 
expect that our parameter estimation for a real dataset 
should be able to determine real /3 p and (3 S with reason- 
able accuracy. 

4. RESULTS 

We display the results of our fitting procedure in Figure 
[8j which shows a map of x 2 distribution as a function of 
j3 p and /3 S . We see that in the broad region near the min- 
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TABLE 1 

DI Herculis parameters 



Parameter 


Primary 


Secondary 


Stellar radius (Rq) 


2.68 


2.48 


Stellar mass (Mq) 


5.15 


4.52 


Von Zeipel's parameter /3 g 


0.1 


0.1 


Effective temperature (K) 


17300 


15400 


v ro t sin (3 (km • s _1 ) 


108 


116 


Cl 


0.35 


0.64 


A(°) 


72 


-84 


Derived parameters 


n 


62 ± 17 


90 < Ps < HO 


ljR (km • s 1 ) 


112 


124 


S 


0.040 


0.046 



imum the value of x 2 is almost constant, which precludes 
us from deriving accurate values of the spin angles from 
the eclipse analysis alone. While the angle (3 P for the 
primary is constrained to lie in the range 30° — 70°, the 
eclipse photometry alone does not set a reasonable limit 
on (3 S : we can only say that it should lie within 50° — 140° 
interval. This difference is caused by the different noise 
levels for primary and secondary eclipses: 4 out of 5 sec- 
ondary eclipses used have dj > 0.01 mag, while 3 out 
of 4 primary eclipses have Uj < 0.01 mag (one primary 
eclipse has aj — 0.004 mag). As we demonstrated in pre- 
vious section, large values of (jj significantly deteriorate 
the performance of our parameter estimation procedure 
(see Figure [7]), which is apparently the case for secondary 
eclipses, during which the lightcurve is most sensitive to 

To obtain a better measurement of these angles we ap- 
ply two additional constraints. One of them uses the ob- 
served precession rate Co b s = 1°. 24 ± 0°. 18/100 yr (Mar- 
tynov & Khaliullin 1980). The apsidal precession rate 
Co depends on /3 PiS since it contains a contribution due 
to the rotation-induced stellar quadrupole, see equation 
(|B18|) , while the latter depends on these angles according 
to equation ([7]). We show the constraint on the apsidal 
precession rate (corresponding to la deviation) by yellow 
curve in Figure [5} w here the analytical estimate of Co is 
obtained using (|B18|) . 

To obtain approximate values and the error bars of 
the angles f3 p and (3 S based on the eclipse fitting and 
the measurement cj, we first construct the photomet- 
ric probability distribution of these angles using the % 2 
map from Figure [5] We then convolve it with the distri- 
bution of P p and p s (assumed to be a two-dimensional 
Gaussian) based on the Co measurement of Albrecht et 
al. (2009). The map of the resultant probability den- 
sity distribution is shown in Figure [51 From this map 
we find /3 p = 62° ±17° and /3 S = 90° ± 20°, where the 
errors correspond to 1-cr uncertainty. Comparing with 
Figure [5] we see that (3 S is constrained essentially purely 
by the Co measurement, with photometric data not be- 
ing useful. At the same time, for the primary angle f3 p 
the photometric data do result in a meaningful measure- 
ment, reducing [3 p from the value suggested by Co alone 
and lowering error considerably. 

Another constraint on spin orientation is based on the 
evolution of v ro t sin /? for both components (see £13.3 1 and 
is illustrated by the white dashed line in Figures [8] & [9j 
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Fig. 8. — x 2 distribution over all eclipses. The black dot shows 
the global minimum of \ 2 distribution within the considered range 
of P and /3 S . The yellow ellipse shows the constrain coming from 
the precession rate corresponding to l a level, where the analytical 
estimate of ill is obtained using l lB18t . Evolution of v ro t sin /3 at 
present time additionally constrains f} p < 90°, /3 S > 90° (repre- 
sented by white dashed lines). 




Fig. 9. — Probability density distribution for /3 P , f} 3 obtained by 
combining the eclipse fitting and the constraint on the precession 
rate. Yellow thick contour gives lcr level. The best fit values of 
spin angles that we derive from this map are /3 P = 62° ± 17° and 
fis = 90° ± 20° (yellow thin contour). The white dashe d cu rve 
additionally shows the v ro t sin/3 evolution constraint, see A3. 31 

This constraint is most important for the spin orienta- 
tion of the secondary as it excludes /3 S < 90° from the 
consideration. As a result, we come up with a refined 
measurement of /3 S = 100° ± 10°. 

Figure [5] shows model lightcurves for these best fit 
values of j3 p and (3 S for all eclipses used in this work. 
One can clearly see the existence of some features in 
the lightcurves that remain unfit by our gravity-darkened 
model, especially at the midpoint of some eclipses. These 



140 
120 
| 100 



20 



A 


— primary 




(a) 




ft 


= 62' 




o 


- secondary 


= 100" 




a 




1 













1950 1960 1970 1980 1990 2000 2010 



140 

120 

I 100 

I 80 
> 

> 40 
20 



A 


- primary 






(b) 






Pp 


= 55° 




o 


- secondary 


A = 


-- 120" 










4 / 















1950 1960 1970 1980 1990 2000 2010 



Fig. 10. — Evolution of Vj sin 0j for both stars for derived angles 
P P = 62°, B = 100° (a) and for another set P = 55°, /3 B = 120° 
(b). In both cases X p = 70° and A s = —76°. Solid and dashed 
curves correspond to theoretical curves for primary and secondary 
stars. Circles and triangles with errorbars (when available) rep- 
resent the measurements for secondary and primary, correspond- 
ingly, taken from Albrecht et al. (2009). It shows that the derived 
/3 S from photometrical analysis does not provide the best fit for 
V s sin (3 S evolution. 

are likely artefacts of the measurements using different 
instruments and at different locations. 

5. DISCUSSION AND CONCLUSIONS 

In this work we developed a method for determining 
full three-dimensional spin-orbit geometry of an eclipsing 
binary system with rapidly rotating components. The 
idea behind this method lies in using the gravity darken- 
ing effect and its influence on the properties of the pho- 
tometric lightcurve of the system. Using this method, 
coupled with two additional constraints — the value of 
the apsidal precession rate of the system and the evo- 
lution of spectroscopically determined projection of the 
stellar rotation speed — we were able to provide a rea- 
sonable measurement of the projections of stellar spins 
onto our line of sight in the eclipsing binary DI Her. 

A very similar technique based on the gravity dark- 
ening effect has already been employed to infer the spin- 
orbit orientation in a planetary system KOI-I3.0I, which 
contains a rapidly rotating (v rot sin/3 = 65 ± 10 km s _1 ) 
intermediate mass star (Szabo et al. 2011; Barnes et al. 
2011). This measurement used several eclipses (transits) 
obtained over a short time span, as opposed to our pro- 
cedure that uses data spread over a long time interval. 
In the case of Barnes et al. (2011) the exquisite photo- 
metric accuracy of Kepler allowed derivation of a rather 



tight constraint on the spin orientation of the host star 
in KOI-13.01 system, something that we cannot accom- 
plish with our low-quality multi-epoch photometry. The 
same kind of photometric accuracy (~ 10~ 4 mag) would 
provide us with a much better constraint on the DI Her 
orientation angles even with single epoch data, see £12.41 

Modeling the gravity darkening-modified eclipse 
lightcurves in binary stars is not an easy task because 
each component of the binary covers large portion of the 
disk of another. This requires integration of the inten- 
sity distribution over a large fraction of the stellar sur- 
face, which naturally gives rise to degeneracies between 
different parameters of the system, making it difficult to 
determine stellar spin orientation. On the other hand, 
in the case of planetary transits planet covers only a 
small fraction of the stellar surface so that the eclipse 
lightcurve can be directly related to the one-dimensional 
run of stellar surface temperature asymmetries. The lat- 
ter can be much more easily modeled via the gravity 
darkening effect to infer the system orientation. 

As of now there is no good explanation for the strong 
spin-orbit misalignment of DI Her (Albrecht et al. 2010). 
It could be primordial, resulting from an interaction be- 
tween the stars and the disk from which they formed, 
which would require some yet unknown mechanism to 
get accomplished. Alternatively, the system may contain 
a third body in a wider orbit as suggested by Kozyreva 
& Bagaev (2009) based on timing of eclipses over a long 
time span. If the orbit of that body is highly inclined 
with respect to the orbit of the inner two stars then 
the Lidov-Kozai mechanism (Lidov 1962; Kozai 1962) 
may be invoked to explain the spin-orbit misalignment 
in DI Her. This idea clearly requires further investiga- 
tion, but we will mention that Lidov-Kozai cycles with 
tidal dissipation are often considered responsible (Fab- 
rycky & Tremaine 2007) for the spin-orbit misalignments 
inferred in many extrasolar planetary systems (Albrecht 
et al. 2012). 

It is interesting that recently started BANANA 
project (Albrecht et al. 2011) focusing on the Rossiter- 
McLaughlin measurements in binaries containing rapidly 
spinning stars reported close spin-orbit alignment for the 
primary star in the NY Cep system, which is similar to DI 
Her in many respects. On the other hand, there are other 
eclipsing binary systems such as AS Camelopardalis, 
which exhibit anomalously slow apsidal precession, simi- 
lar to DI Her (Pavlovski et al. 2011). If it is the spin-orbit 
misalignment that is causing the anomalous precession 
in AS Camelopardalis then the photometric method de- 
veloped in this work and Barnes (2009) may be used to 
constrain the system orientation (although the measured 
projected rate in this system is not very high, 15 km 
s _1 for the primary component). 

Authors are grateful to the referee, Jason Barnes 
for valuable suggestions that helped us improve the 
manuscript, Valentina Kozyreva for providing us with 
observational data, and Ed Turner for useful discussions. 
AAP thanks Princeton University for hospitality during 
the time when part of this work was performed. Finan- 
cial support of this research is provided by the Dinasty 
fellowship for AAP and by the Sloan Foundation and 
NASA via grant NNX08AH87G for RRR. 
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APPENDIX 

GEOMETRY OF THE SKY-PROJECTED STELLAR DISC 



In this section we will derive the relation between the coordinates of a given point on the stellar surface in the 
observer and symmetry frames. Normal £ to the stellar surface at a point (x, y, z) in the symmetry frame is (|£| ^ 1) 



T} 2 ' f] 2 



(Al) 



Because of the rotation-induced distortion the sky-plane projected stellar shape is not circular. The last visible points 
on the stellar surface for the observer are given by the following equation 



n-£ = 

where n is a unit vector towards the observer given by equation It results in the following equations: 



tan f3 



-x, 



1 



tan 2 /? 



y 



= R 



pol 



(A2) 



(A3) 



Equations (|A3j) give us the " critical line" of the last visible points on the stellar surface (where £ lies in the sky plane) 
in the symmetry system. It defines the shape of the visible stellar disk. To obtain the coordinates in observer frame 
the corresponding coordinate transformation should be made: 



x = cos /3xq — sin/3Z, z = sin/3xo + cos 

where we defined 

x = cos XX + sin XY, y = — sin AA" + cos XY. 

So the equation for the critical line written in the observer frame is 

x (r]~ 2 — 1) tan/3 



1 + ij- 2 tan 2 (3 



Xq tanp, 



r 2 

•' o 



2 

Vo 



T] Rp 0l 



cos 2 (3 + 1] 2 sin j3 

where p is defined by equation (IA6|) . In coordinates (xq, yo, Z) the stellar surface is described by 

Uo + (cos(/3)x - sin(/3)Z) 2 



+ (sm(P)x + cos(P)Zy 



r>2 



from which one can find Z in terms of X and Y . The solution to the resulting quadratic is 

~ x -fdet 



Z = x^ tanp + 



V 



cos 2 /? + r\ 2 sin (3 



where the expression for the determinant det is 



det = 



+ (yl - V 2 R 2 pol ) cos 2 (3 + 



sin 2 /3 



(A4) 
(A5) 

(A6) 
(A7) 

(A8) 

(A9) 

(A10) 



It can be easily checked that the condition det = coincides with the equation for the critical line, so in the region 
interior to this line det > 0. We choose the positive root of the determinant (the negative root represents the invisible 
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side of the star as seen from Earth, see Barnes (2009). So if (X, Y, 0) is the point on the sky-projected stellar disc then 
R(X, Y) = (X,Y, Z(X,Y)) represents the full three-dimensional coordinates of the point on the stellar surface right 
above (X, Y, 0) in the observer frame. 

In general, to find coordinates (x, y, z) of a point on the stellar surface in the symmetr y fra me corresponding to 
a poi nt (X , Y) projected onto the sky plane one has to compute xq and y$ via equation (|A5|) . determine Z using 
(|A9j) - (|A10j) . obtain x and z from equation (|A4j) . and finally determine y from equation ([5]). 

TIME EVOLUTION 

Here we derive equations that describe the evolution of the binary orbit and spin orientation of its components. 
We denote L = LI, Sj = SjSj, j = p,s the orbital angular momentum and spin angular momenta of the two stars 
respectively, |1| = |s 3 -| = 1. Here 



L = /iO^-a 2 -\/l — e 2 , Sj = IjU>j = r/jMjR^ujj, (Bl) 

where Mj, Rj, Wj are the masses, radii and spin angular frequencies of the two stars, rjj are constants determining 
their moments of inertia, /i = M p M s /(M p + M s ) is the reduced mass, a, e, and fix = [G(M p + M s ) / a 3 ] 1//2 are the 
semi-major axis, eccentricity and mean orbital frequency. We introduce unit vector n from the system's barycenter to 
the observer and direct a Cartesian coordinate system in the directions n, m, k where 

k=^^, m = kxn= — [l-n(l-n)l, (B2) 
sin i sin i 

where i is the inclination of the system, cosi = (1 • n), see Figure [T] 

Orientation of Sj is conventionally given by the three angles ctj, (3j, between each of Sj and vectors 1, n, and k 
correspondingly, i.e. 

cosaj = (sj ■ 1), cos/3j = (sj ■ n), cos7j = (sj ■ k). (B3) 

Then it is easy to show that 

cos a,- — cosicos/3,- , 

s, = cos PjU H + cos 7,k. (B4j 

sini 

Observationally, it is also convenient to introduce angle between the projection s±j — [sj — n(sj • n)]/ sin f3j of s^- 
onto the sky plane and the projection m of vector 1 onto the same plane: cosAj = (s±.j ■ m). This is the angle which 
is measured by the Rossiter-McLaughlin effect. One can easily show that these four angles are related via 

cos ctj = cos i cos fij + sin i sin j3j cos Xj , cos 7 = sin A sin /?. (B5) 

Thus, knowing Aj, (3j and i one can immediately obtain ctj from these expressions. 

Orientation of the orbital ellipse in the plane of the orbit is given by the eccentricity vector E = ee which points 
from the main focus to the pericenter. Given that direction defined by vector k is the direction of the line of nodes we 
identify the angle between k and e as the longitude of the periastron ui and write 

e = cos wk + sin w(l x k). (B6) 

Stellar asphericity due to rotation and tides as well as relativistic effects lead to evolution of 1, e, and Sj described 
by the following equations (Barker & O'Connell 1975): 

1 = fij, x 1, e = fij, X e, Sj = ils,j x Sj, (B7) 

where 



n L = ln E + n T ,j 1 + n < 



1 — 5 cos ctj , 
cosa,s,- H M 



(B8) 



L Q,s 

Sls,j = ^Gjl + &p,j i s j — 3cosajl) ~ f2p,j {sj ~ 3cosajl) . (B9) 

Here different contributions to precession rates are denoted as follows (Barker & O'Connell 1975; Claret et al 2010): 
orbital Einstein precession 

3GCl K (M s + M P ) 
c 2 a(l — e 2 ) 

orbital precession caused by stellar quadrupole due to tidal distortions (&2j are introduced below) 



Me = 0JGR = 2 2 , P ~ « 2.35 o /100 yr, (B10) 



M r {RjY 8 + 12e 2 



^ = 15 ^^7j W r ^ ij (Bn) 

tt T , p ~ 0.69°/100 yr, Q T ,s ~ 0.637100 yr, 
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orbital precession due to rotation-induced stellar quadrupole 



3 G(M p + M s )J 2d _ M p + M, 



2 ft*ra 5 (l-e 2 ) 2 



2.1 o /100 yr 

(sin/3 p ) 2 ' 



Mj n K (l-e 2 ) 2 V a 
2.2 o /100 yr _ x 

(sin/3 fl )2 S ' 



(B12) 



geodetic spin precession 



, 3 ^ r, 



2c 2 a(l - e 2 ) 
G , P fa 0.617100 yr, G , S w 0.697100 yr, 

and the spin precession caused by rotation-induced stellar oblateness 



(B13) 



GM p M s J 2 ,i 



k 2> j M r ujj 



n 



p,p 



2/ j w 3 a 3 (l - e 2 ) 3 / 2 3% Mj (1 - e 2 ) 3 / 
136.67100 yr 



(B14) 



168.67l00yr _ a 

5 i 



sin/3 p sin/3 s 

In equations (IB11[) . (|B12I) . (|B14I) fc 2j - is the apsidal motion constant related to stellar rotation-induced quadrupole 
moment constant J 2 .j via 

1/2 



3 J 2 j I Ldf, 



k2 ' 3 ~ 2 m ( u < 



(B15) 



with u)h.j being the breakup angular frequency. Numerical estimates assume k 2 j ps 0.008 (Claret et al. 2010) and 
the moment of inertia co nsta nt rjj = 0.063 (Claret & Gimenez 1989). Given that &G,j "C &p,j we dropped geodetic 
contribution in equation (|B9j) . 
Differentiating relation cosi = (1-n) with respect to time and using equations (|B2|) . (|B3|) . (|B7[) . and (IB8|) one obtains 



In 



i = — : — = SIl • k = fin cos a, cos7,-. 

sin r ^ ' J J J 



(B16) 



Because of precession of 1 vectors k and m vary in time. Differentiating equations (|B2[) with respect to time and using 
(|B3|) . (1B7[) . and (|B8[) their evolution is governed by equations 



. 2 

sm 1 



57q j cos aj (cos i cos ay — cos f3j ) 



3=P,s 



sin 2 1 



2_, ^Iqj cos a j (cos (3 j — cos i cos a j) , (B17) 



Next, we differ enti a te co s oj = ( e ■ k ) as w ell as each of the relations (|B3|) with respect to time and transform them 
using equations (IB2|) - (IB4|1 . (|B6[) - (|B9[) . (IB16|) and (IB17|) . As a result we arrive at the following expressions: 



u - 



03 



(e • k) + (e • k) 



smw 



= =n E + J2 \ q t,j 



sin 2 i 



cos a.j (cos a j — cos i cos /3j ) + sin i 



, .1 — 5 cos 2 <x 3 - 



(1-Si) + (1-Bj) 



=-n £ 



sm otj 

cosa r 
sin i sin a. 



[cos 7 r (cos i cos ctj — cos (3j) — cos 7j (cos i cos a r — cos j3 r )] , j ^ 



(kj ■ n) sin i cos 7,- cos a, 
— = 3s2pj- 



sin j3j 

(k ■ Si ) + (k ■ Sj) 

sin 7,; 
cos ay — cos i cos j3j 



sin (3j 



+ 



sm i sm 7o 

'J r—s,p 

3£1 p j cos ay (cos i cos a 3 - — cos f3j ) 
sin i sin 7^ 



flQ >r cos a r (cos i cos a r — cos (3 r ) 



(B18) 

(B19) 
(B20) 



(B21) 
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Equations (IB16[) . (|B18I) - (|B21|) constitute a closed system of 8 evolution equations for 8 unknown angles — i, u, ctj, 
7j, j — p, s — fully determining the orbital orientation of the binary and spin orientation of each star. 
One can check the validity of these expressions by using the fact that the total angular momentum of the system 
J = LI + SpSp + S s s s is conserved. Equation (|B18|) agrees with the analogous expression in Shakura (1985). 



